################################
#MAPS HABSBURG EMPIRE CHOROPLETH#
################################
#This to obtain:
#Figure_2a.jpg
#Figure_2b.jpg
#Figure_2c.jpg
#Figure_2d.jpg

rm(list = ls())
library("ggplot2")
library("dplyr")
library("ggrepel")
library("raster")
library("RColorBrewer")
library("lemon")
library("sp")
library("maptools")
library("ggsn")
library("ggrepel")
library("sf")
library("ggplot2")
library("stringr")
library("readxl")

#setwd("C:/Users/bogdanp/")
setwd("/Users/bgpopescu/")

#Reading Polygons
HABS <- st_read(dsn="./Dropbox/Legacies_Project/Analysis/data/data.gdb",
                    layer="all_hasburg_GCS")
sf::sf_use_s2(FALSE)
HABS<-st_simplify(HABS,  dTolerance = 0.005)

all_with_data_croatia <- st_read(dsn="./Dropbox/Legacies_Project/Analysis/data/data.gdb",
                layer="all_with_data_croatia_GCS")
all_with_data_croatia<-st_simplify(all_with_data_croatia,  dTolerance = 0.005)

mil_col_short_polygon <- st_read(dsn="./Dropbox/Legacies_Project/Analysis/data/data.gdb",
                  layer="krajna_line_GCS")
mil_col_short_polygon<-st_simplify(mil_col_short_polygon,  dTolerance = 0.005)

mil_col_polygon <- st_read(dsn="./Dropbox/Legacies_Project/Analysis/data/data.gdb",
                              layer="krajna_polygon_d_GCS")
mil_col_polygon<-subset(mil_col_polygon, select = c(name_regiment))
mil_col_polygon<-st_simplify(mil_col_polygon, dTolerance = 0.005)

all_coun<- st_read(dsn="./Dropbox/Legacies_Project/Analysis/data/data.gdb",
                   layer="all_countries_GCS")
all_coun<-st_simplify(all_coun, dTolerance = 0.005)

all_coasts<- st_read(dsn="./Dropbox/Legacies_Project/Analysis/data/data.gdb",
                   layer="only_coasts_GCS")
all_coasts<-st_simplify(all_coasts, dTolerance = 0.005)

#Reading Rasters
euro_raster <- raster("./Dropbox/Legacies_Project/Analysis/data/euro_raster.tif")

#Saving coordinates for the Habsburg polygon
HABS$y_lat<-sf::st_coordinates(st_centroid(HABS))[,2]
HABS$x_lon<-sf::st_coordinates(st_centroid(HABS))[,1]

max_lat<-max(HABS$y_lat) + 4
min_lat<-min(HABS$y_lat) - 4

max_lon<-max(HABS$x_lon) + 4
min_lon<-min(HABS$x_lon) - 4

e <- extent(min_lat, max_lat, min_lon, max_lon)
rc <- crop(euro_raster, e)	
res(rc)

#Lower resolution raster
rc.aggregate <- aggregate(rc, fact=5)
res(rc.aggregate)

#Get rid of really high mountains
rc <- reclassify(rc.aggregate, c(10000,Inf,NA))
rm(rc.aggregate)
rna <- reclassify(rc, cbind(NA, 0))
rm(rc)
#Convert rasters TO dataframes for plotting with ggplot
hdf <- rasterToPoints(rna); hdf <- data.frame(hdf)
colnames(hdf) <- c("X","Y","Elev")

#   Create vectors for colour breaks
b.hs <- seq(min(hdf$Elev),max(hdf$Elev),length.out=100)

# create slope and hillshade
slope = terrain(rna, opt='slope')
aspect = terrain(rna, opt='aspect')
hill = hillShade(slope, aspect, 40, 270)

#######################
#Creating color breaks#
#######################

no_classes <- 6
labels <- c()
HABS$pct_military<-(HABS$military/HABS$population)*100

obj <- c(0, 0.3, 0.5, 0.6, 1, 5, 10.4575)
quantiles<-obj

# here I define custom labels (the default ones would be ugly)
labels <- c()
for(idx in 1:length(quantiles)){
  labels <- c(labels, paste0(round(quantiles[idx], 2), 
                             " - ", 
                             round(quantiles[idx + 1], 2)))
}
# I need to remove the last label 
# because that would be something like "66.62 - NA"
labels <- labels[1:length(labels)-1]

# here I actually create a new 
# variable on the dataset with the quantiles
HABS$brks <- cut(HABS$pct_military, 
                    breaks = quantiles, 
                    labels = labels, 
                    include.lowest = T)

HABS$latlon<-st_centroid(st_geometry(HABS))
HABS[c('lat', 'lon')]<-str_split_fixed(HABS$latlon, ", ", 2)
HABS$lat<-as.numeric(str_replace(HABS$lat, "c\\(", ""))
HABS$lon<-as.numeric(str_replace(HABS$lon, "\\)", ""))

data_1851<-read_excel("./Dropbox/Legacies_Project/Analysis/data/census_1851_general.xlsx", sheet=1, col_names = TRUE, skip = 0)
data_1851<-subset(data_1851, select = c(name_map_english, houses, familien, persons_per_family))
HABS<-left_join(HABS, data_1851, by = c("name_english"="name_map_english"))

lbl <- subset(HABS, select = c(name_english, lat, lon))
my.cols <- brewer.pal(6, "Greys")

map<-ggplot() + 
  scale_alpha(name = "Altitude in m", guide="none")  + 
  geom_sf(data = HABS, colour = "grey60", aes(fill=brks))+
  geom_sf(data = mil_col_short_polygon, colour = "blue", fill = NA, linewidth = 1)+
  geom_sf(data = all_coasts, colour = "black", fill = NA, linewidth = 0.5)+
  coord_sf(ylim=c(40.815009,53), xlim=c(8.556495, 27))+
  theme_bw()+
    theme(plot.title = element_text(hjust = 0.5))+
    scale_fill_manual("", values = my.cols, guide = "legend") +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.title.y=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank(),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.title=element_text(size=14),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12),
        aspect.ratio = 1)+
  geom_label_repel(data=lbl, aes(label=lbl$name_english, x = lat, y = lon,
                                 fontface = 'bold'), alpha = 0.7)
map

map2<-map+ggsn::scalebar(x.min = 10, x.max = 14,
           #Above you mention how long the scalebar should be
           y.min = 52, y.max = 53,
           #Above you mention how thick the scalebar should be
           dist = 200, dist_unit = "km",
           transform = T, model = "WGS84",
           location = "topright",
           st.size = 4,
           #Above you have the font size of the numbers below the scalebar
           st.dist =0.2,
           #Above you have the distance between the bar and the text, as a proportion of the y axis.
           height=0.2)
# #Above you have anumber between 0 and 1 to indicate the height 
# #of the scale bar, as a proportion of the y axis
map3<-reposition_legend(map2, 'bottom left')

ggsave(map3,file="./Dropbox/Legacies_Project/Paper/figures/Figure_2a.jpg", 
       height=19.42, width=20.2, units = "cm", dpi=300)



#Persons per Family
no_classes <- 6
labels <- c()

obj <- c(2, 4, 5, 6, 7, 8, 9)
quantiles<-obj

# here I define custom labels (the default ones would be ugly)
labels <- c()
for(idx in 1:length(quantiles)){
  labels <- c(labels, paste0(round(quantiles[idx], 2), 
                             " - ", 
                             round(quantiles[idx + 1], 2)))
}
# I need to remove the last label 
# because that would be something like "66.62 - NA"
labels <- labels[1:length(labels)-1]

# here I actually create a new 
# variable on the dataset with the quantiles
HABS$brks <- cut(HABS$persons_per_family, 
                 breaks = quantiles, 
                 labels = labels, 
                 include.lowest = T)

map<-ggplot() + 
  scale_alpha(name = "Altitude in m", guide="none")  + 
  geom_sf(data = HABS, colour = "grey60", aes(fill=brks))+
  geom_sf(data = mil_col_short_polygon, colour = "blue", fill = NA, linewidth = 1)+
  geom_sf(data = all_coasts, colour = "black", fill = NA, linewidth = 0.5)+
  coord_sf(ylim=c(40.815009,53), xlim=c(8.556495, 27))+
  theme_bw()+
  theme(plot.title = element_text(hjust = 0.5))+
  scale_fill_manual("", values = my.cols, guide = "legend") +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.title.y=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank(),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.title=element_text(size=14),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12),
        aspect.ratio = 1)+
  geom_label_repel(data=lbl, aes(label=lbl$name_english, x = lat, y = lon,
                                 fontface = 'bold'), alpha = 0.7)
map

map2<-map+ggsn::scalebar(x.min = 10, x.max = 14,
                         #Above you mention how long the scalebar should be
                         y.min = 52, y.max = 53,
                         #Above you mention how thick the scalebar should be
                         dist = 200, dist_unit = "km",
                         transform = T, model = "WGS84",
                         location = "topright",
                         st.size = 4,
                         #Above you have the font size of the numbers below the scalebar
                         st.dist =0.2,
                         #Above you have the distance between the bar and the text, as a proportion of the y axis.
                         height=0.2)
# #Above you have anumber between 0 and 1 to indicate the height 
# #of the scale bar, as a proportion of the y axis
map3<-reposition_legend(map2, 'bottom left')

ggsave(map3,file="./Dropbox/Legacies_Project/Paper/figures/Figure_2c.jpg", 
       height=19.42, width=20.2, units = "cm", dpi=300)



######################
#MAP Specific Croatia#
######################

all_with_data_croatia$y_lat<-sf::st_coordinates(st_centroid(all_with_data_croatia))[,2]
all_with_data_croatia$x_lon<-sf::st_coordinates(st_centroid(all_with_data_croatia))[,1]

max_lat<-max(all_with_data_croatia$y_lat) + 4
min_lat<-min(all_with_data_croatia$y_lat) - 4

max_lon<-max(all_with_data_croatia$x_lon) + 4
min_lon<-min(all_with_data_croatia$x_lon) - 4


#######################
#Creating color breaks#
#######################

no_classes <- 6
labels <- c()

names(all_with_data_croatia)
obj <- c(0, 0.3, 0.5, 0.6, 1, 5, 10.4575)
quantiles<-obj

# here I define custom labels (the default ones would be ugly)
labels <- c()
for(idx in 1:length(quantiles)){
  labels <- c(labels, paste0(round(quantiles[idx], 2), 
                             " - ", 
                             round(quantiles[idx + 1], 2)))
}
# I need to remove the last label 
# because that would be something like "66.62 - NA"
labels <- labels[1:length(labels)-1]

# here I actually create a new 
# variable on the dataset with the quantiles
all_with_data_croatia$brks <- cut(all_with_data_croatia$pct_military, 
                    breaks = quantiles, 
                    labels = labels, 
                    include.lowest = T)

sf_use_s2(FALSE)


all_with_data_croatia$latlon<-st_centroid(st_geometry(all_with_data_croatia))
all_with_data_croatia[c('lat', 'lon')]<-str_split_fixed(all_with_data_croatia$latlon, ", ", 2)
all_with_data_croatia$lat<-as.numeric(str_replace(all_with_data_croatia$lat, "c\\(", ""))
all_with_data_croatia$lon<-as.numeric(str_replace(all_with_data_croatia$lon, "\\)", ""))
lbl <- subset(all_with_data_croatia, select = c(name_merge, lat, lon))

my.cols <- brewer.pal(6, "Greys")
all_with_data_croatia$brks


map2<-ggplot() + 
  scale_alpha(name = "Altitude in m", guide="none")  + 
  geom_sf(data = all_with_data_croatia, colour = "grey60", aes(fill=brks))+
  geom_sf(data = mil_col_polygon, colour = "blue", fill = NA, linewidth = 1.5)+
  geom_sf(data = all_coun, colour = "black", fill = NA, linewidth = 0.5)+
  #This is for Croatia
  coord_sf(ylim=c(42.4,46.669853), xlim=c(13.556495, 19.447201))+
  theme_bw()+
  scale_fill_manual("", values = my.cols, guide = "legend") +
  geom_label_repel(data=lbl, aes(label=name_merge, x = lat, y = lon,
                                 fontface = 'bold'), alpha = 0.7)+
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.title.y=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank(),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.title=element_text(size=14),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12),
        aspect.ratio = 1)+
  ggsn::scalebar(x.min = 18.5, x.max = 19.1,
                 #Above you mention how long the scalebar should be
                 y.min = 46.4, y.max = 46.7,
                 #Above you mention how thick the scalebar should be
                 dist = 50, dist_unit = "km",
                 transform = T, model = "WGS84",
                 location = "topright",
                 st.size = 4,
                 #Above you have the font size of the numbers below the scalebar
                 st.dist =0.2,
                 #Above you have the distance between the bar and the text, as a proportion of the y axis.
                 height=0.2)

map2<-reposition_legend(map2, 'bottom left')

ggsave(map2,file="./Dropbox/Legacies_Project/Paper/figures/Figure_2b.jpg", 
       height=19.42, width=20.2, units = "cm", dpi=300)


excel_1851 <- read_excel(".//Dropbox/Legacies_Project/Analysis/data/census_1851.xlsx", sheet=1, col_names = TRUE, skip = 0)
all_with_data_croatia_1851 <- st_read(dsn="./Dropbox/Legacies_Project/Analysis/data/data.gdb",
                                      layer="all_with_data_croatia_1851_GCS")
all_with_data_croatia_1851<-st_simplify(all_with_data_croatia_1851,  dTolerance = 0.005)
data_1851<-left_join(all_with_data_croatia_1851, excel_1851, by = c("merge_name"="name_merge"))


#######################
#Creating color breaks#
#######################

no_classes <- 4
labels <- c()

my.cols <- brewer.pal(4, "Greys")
quantiles<-quantile(data_1851$persons_per_family, na.rm=T)

# here I define custom labels (the default ones would be ugly)
labels <- c()
for(idx in 1:length(quantiles)){
  labels <- c(labels, paste0(round(quantiles[idx], 2), 
                             " - ", 
                             round(quantiles[idx + 1], 2)))
}
# I need to remove the last label 
# because that would be something like "66.62 - NA"
labels <- labels[1:length(labels)-1]
labels

# here I actually create a new 
# variable on the dataset with the quantiles
data_1851$brks <- cut(data_1851$persons_per_family, 
                      breaks = quantiles, 
                      labels = labels, 
                      include.lowest = T)

data_1851_b<-subset(data_1851, !is.na(data_1851$persons_per_family))

data_1851_b$lon<-st_coordinates(st_centroid(data_1851_b))[,1]
data_1851_b$lat<-st_coordinates(st_centroid(data_1851_b))[,2]
data_1851$name_merge


all_coun <- st_read(dsn="./Dropbox/Legacies_project/Data/Historical Data - Habsburg Empire/final.gdb",
                    layer="all_countries_GCS")
sf::sf_use_s2(FALSE)
all_coun<-st_simplify(all_coun,  dTolerance = 0.005)


lbl <- subset(data_1851_b, select = c(name_merge, lat, lon))
#Removing some labels for specific districts
lbl$name_merge[lbl$name_merge=="Petrinja"]<-""
lbl$name_merge[lbl$name_merge=="Kreuz"]<-""
lbl$name_merge[lbl$name_merge=="Belovar"]<-""
lbl$name_merge[lbl$name_merge=="Ivanic"]<-""
lbl$name_merge[lbl$name_merge=="Carlopago"]<-""
lbl$name_merge[lbl$name_merge=="Zengg"]<-""
lbl$name_merge[lbl$name_merge=="Brod"]<-""


mil_col_polygon <- st_read(dsn="./Dropbox/Legacies_Project/Data/Maps/empires.gdb",
                              layer="krajna_polygon_d_GCS")
mil_col_polygon<-subset(mil_col_polygon, select = c(name_regiment))
mil_col_polygon<-st_simplify(mil_col_polygon, dTolerance = 0.005)


map<-ggplot() + 
  scale_alpha(name = "Altitude in m", guide="none")  + 
  geom_sf(data = all_coun,fill = NA, colour = "grey60")+
  geom_sf(data = data_1851_b, colour = "grey60", aes(fill=brks))+
  geom_sf(data = mil_col_polygon, colour = "blue", fill = NA, linewidth = 1.5)+
  #This is for Croatia
  coord_sf(ylim=c(42.4,46.669853), xlim=c(13.556495, 19.447201))+
  theme_bw()+
  labs(x = "Longitude", y="Latitude")+
  #ggtitle("Family Size, 1851")+
  theme(plot.title = element_text(hjust = 0.5))+
  scale_fill_manual("", values = my.cols, guide = "legend") +
  geom_label_repel(data=lbl, aes(label=lbl$name_merge, x = lon, y = lat,
                                 fontface = 'bold'), alpha = 0.7)+
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.title.y=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank(),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.title=element_text(size=14),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))+
  ggsn::scalebar(x.min = 18.5, x.max = 19.1,
                 #Above you mention how long the scalebar should be
                 y.min = 46.4, y.max = 46.7,
                 #Above you mention how thick the scalebar should be
                 dist = 50, dist_unit = "km",
                 transform = T, model = "WGS84",
                 location = "topright",
                 st.size = 4,
                 #Above you have the font size of the numbers below the scalebar
                 st.dist =0.2,
                 #Above you have the distance between the bar and the text, as a proportion of the y axis.
                 height=0.2)

map2<-reposition_legend(map, 'bottom left')

ggsave(map2,file="./Dropbox/Legacies_Project/Paper/figures/Figure_2d.jpg", 
       height=19.42, width=20.2, units = "cm", dpi=300)
